clear; clc; close all;


mu = -0.6;
k_fixed = 1.78;
x0_fixed = 0.8;
q0_range_top = linspace(-15,15,2000);  % 上图 q0 范围
transient_top = 300;
record_top = 100;

% 用于下图：在 (x0, q0) 平面上扫描初始条件
Nx = 200;
Nq = 200;
x0_range = linspace(-2,2,Nx);
q0_range = linspace(-15,15,Nq);
transient_bottom = 100;
total_iter_bottom = 1000;  % 包含 transient_bottom


%% 上图：固定 x0 = 0.8，扫描 q0，绘制 q_final vs. q0
q_final_all = zeros(1, length(q0_range_top));
for i = 1:length(q0_range_top)
    q0 = q0_range_top(i);
    % 初始条件 (x0_fixed, q0)
    x = x0_fixed;
    q = q0;
    % 舍去 transient_top 步
    for t = 1:transient_top
        [x, q] = MTM(x, q, mu, k_fixed);
    end
    % 记录 record_top 步后的 q 平均值（或最后一个 q）
    q_vals = zeros(1, record_top);
    q_vals(1) = q;
    for t = 2:record_top
        [x, q] = MTM(x, q, mu, k_fixed);
        q_vals(t) = q;
    end
    % 这里简单取最后一个 q 作为代表
    q_final_all(i) = q_vals(end);
end


%% 下图：局部吸引域图，扫描 (x0,q0) 平面
basin_map = zeros(Nq, Nx);  % 行：q0, 列：x0
unbounded_threshold = 1e3;
for ix = 1:Nx
    for iq = 1:Nq
        x = x0_range(ix);
        q = q0_range(iq);
        diverge = false;
        for t = 1:total_iter_bottom
            [x, q] = MTM(x, q, mu, k_fixed);
            if abs(x) > unbounded_threshold || abs(q) > unbounded_threshold
                diverge = true;
                break;
            end
        end
        if diverge
            basin_map(ix, iq) = 0; % 0 表示发散（白色）
        else
            % 简单分类：如果最终 x 接近 0，则判为稳定点，否则按 q 分区
            if abs(x) < 1e-4
                basin_map(ix, iq) = -1;  % 稳定点，特殊标记
            else
                % 按 q mod 2π 分类（这里简单取 round(q/(2*pi)) 作为分类）
                m = round(q/(2*pi));
                % 限定 m 在 [-3,2]
                m = max(min(m,2),-3);
                % 将 m 映射到 1~6 （-3→1, -2→2, …, 2→6）
                basin_map(ix, iq) = m + 4;
            end
        end
    end
end


%% 绘图
% 上图：q vs. q0（固定 x0）
subplot(2,1,1);
plot(q0_range_top, q_final_all, 'b.', 'MarkerSize', 10);
xlabel('q_0'); ylabel('q_{final}');
title(['周期性平移 q 效果 (x_0 = ', num2str(x0_fixed), ')']);
grid on;

% 下图：局部吸引域 in (x0, q0) 平面
subplot(2,1,2);
imagesc(x0_range, q0_range, basin_map);
axis xy;  % 保证 y 轴从下到上递增
colormap(jet(8)); % 这里 jet(8) 包括了 -1, 0, 1..6 的色阶，可根据需要调整
colorbar;
xlabel('q_0'); ylabel('x_0');
title(['局部吸引域 (\mu = ', num2str(mu), ', k = ', num2str(k_fixed), ')']);


%% 2D-MTM 映射函数
function [x_next, q_next] = MTM(x, q, mu, k)
    if x < 0.5
        x_next = mu*x + k*cos(q)*x;
    else
        x_next = mu*(1 - x) + k*cos(q)*x;
    end
    q_next = q + x;
end
